JOURNAL OF RESEARCH of the National Bureau of Standards- B. Mathematical Sciences 
Vol. 78B. No. 4, October-December 1974 

Automatic Computing Methods for Special Functions. 
Part II. The Exponential Integral f„(x) 

Irene A. Stegun and Ruth Zucker 

Institute for Basic Standards, National Bureau of Standards, Washington, D.C. 20234 

(September 4, 1974) 

Accurate, automatic, efficient methods for computing the exponential integral E n (x) are detailed and 
implemented in an American National Standard FORTRAN program. The driver program and test 
results are also included. 
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1. Introduction 

The exponential integral itself occurs in many physical problems and often many other inte- 
grals are expressible in terms of the exponential integral. In view of its importance and the dif- 
ficulties accompanying the repeated application of its recurrence relation, we have chosen this inte- 
gral as Part II. (For Part I, see J. of Research NBS, Vol. 74B, July-September 1970, pp. 211-224.) 

Accurate, efficient, automatic computing methods implemented in American National Standard 
FORTRAN, covering the entire range of machine legitimate arguments and/or functional values 
will be supplied. The number of terms in series, the number of convergents in an iterative process, 
etc., are all determined by the computer as a function of word length, arguments, the accuracy 
desired, etc. In cases of error returns, more realistic results will be returned. To further ensure 
correct limiting values of related functions, the proper analytic behavior of the function will always 
be retained. The driver program and test results are also included. 

2. Mathematical Formulas 

Relevant formulas are collected here for completeness and ease of reference. In keeping with 
the convention of the Handbook [1], 1 x here is a real variable. 

A. Definition 
E n (x) = I *t- n e*<dt (n = 0, 1, 2, . . .; x > 0) 



= | u n - 

= x"~ l I t"e 



n - 2 e- x ' u du 

J 



( dt. 



AMS Subject classification: 33-04, 65D20. 
1 Figures in brackets indicate references on page 205. 
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B. Series Expansion 

E n (x)= *' [-Irix + nn)] ~ £' i ( Tu (* > °) 

v ' m=0 v 7 

*(1) =~y, *(n) =-7+ Y - (n>l) 



y (Euler's constant) = .57721 56649 . 



C. Continued Fraction 
1 n 1 n + 1 2 



„ , x /In 1/1 + 12 \ 

^W = ^feiT^-iT-^--J ( * >0) - 

D. Asymptotic Expansion 
E„( x )~ e -\l-^ + n -±±±> _ n(n + l)(n + 2) + 1 

E. Special Values 



£n(0)= V (»> 1) 

Al — 1 



£oU) = — • 

X 



F. Recurrence Relation 



6 X 

E n+ i(x)=ji jiEn(x) 71=1,2, 3,. . .. 



G. Differentiation Formula 

d k 
— E n (x) = (--iyEn- k {x) rc=l,2, 3, . . , 

H. Inequalities 
(* > 0, n= 1, 2, 3, . . .) 

£„(*) < E n+1 (x) < E n (x) 



<e*E„(x) 



x + n x + n — 1 
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r(l-n,x)=x 1 - n E n (x) 
U(n, n, x)=e x x l -"E,,(x). 



I. Related Functions 
Incomplete Gamma Function 

Confluent Hypergeometric Function 

3. Method 



Examination of the recurrence relation indicates that an independent computation must be 
carried out for at least one value of n^\. The function E,,(x) is always positive, for real x has no 
zeros, and is a slowly decreasing function with increasing n and/or x. Repeated application of the 
recurrence relation therefore will yield an increasing round off error even if, with more complicated 
scaling, both forward and backward recurrence relations are used starting at n= [x]. Generally, 
in physical applications the first few orders only are needed. Consequently, we have chosen an 
independent computing method valid for all orders. 

The implementing American National Standard FORTRAN program has been set up in such a 
way as to require a minimal number of changes for varying precision, either single or double pre- 
cision on the same or different computers. The program checks for positive arguments n (=RN) and 
x and integer n. If either is negative, an error indicator is set (IERR = 1) and an impossible value, 
the negative of the maximum machine value (=RINF), is returned for both functional values, 
E„{x) (=ENX) and e J E„(x) (=EXPENX). If n > RMAXI, the maximum integer convertible to a 
floating point number, n is assumed to be an integer. If n ^ RMAXI, an integer test is applied to RN. 
If it fails the test, IERR is set equal to 2, and RINF is returned for both functional values. To assure 
ready portability, additional tests are performed for negative zero and allowance for round off errors 
due either to machine arithmetic or the system being used. 

The special cases are treated independently. When n = 0, to avoid machine difficulty, x is 
tested against the reciprocal of the maximum machine value; if x equals, or is less than that number, 
the maximum value is returned for both functional values. If x is greater, then ENX = e -J 7jt and 
EXPENX = llx. When x = 0, and n = or 1 , ENX = EXPENX - RINF; for n > 1 , ENX = EXPENX 

= 1/(71-1). 

Computation with the series results in a greater round-off error in the neighborhood of x= 1, 
in particular for small values of n. Even if the form of the series were changed from the present alter- 
nating form, there is still a loss due to the logarithmic term. While the continued fraction is con- 
vergent for x > 0, the number of terms increases rapidly as x approaches 0. The integrand is not 
particularly amenable to automatic numerical integration by any of the low order formulas. In 
general (n ^ RMAXI) the most accurate and efficient methods of automatic computation are the 
alternating power series for x ^ l(ULPS) and the "even" form of the continued fraction for x > 1. 

For 1 ^ n ^ RMAXI and < x ^ 1, the power series is used in the following form 

£„(*)= ENX =2 Tm=SVM + TM(M = 0, 1, . . .,£) 
\/ = o 

-{-x) M PTERM 

wheTeTx, = M\tM-n+l = D) = —D forM^-1 

and T M = PTERM [LOG(x) -PSI(n)] forM=n-l. 

PTERM is obtained by recurrence with PTERM(0) = -1 and PTERM(M+1) = 

PTERM (Af) • (-x)/M+l. 
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Rather than compute a previously determined fixed number of terms, it is valid, provided 
M^/i-1, to terminate the series when the relative error, computed as |77lf/SUM| is less than or 
equal to a prescribed tolerance (=TOLER). In the present program, TOLER = 2~ NBM , where NBM 
is the maximum number of binary digits in the mantissa of a floating point number. It may, however, 
be set to the accuracy desired. 

Since the series is an alternating one, SUM may equal zero at some M. However, since the func- 
tion has no zeros for real x 9 the relative error test may be safely bypassed if SUM = and additional 
terms computed. In our manner of evaluating the series, TM may likewise be zero if log(x) = ^P (n). 
(Values of e^ {n) for n— 1(1)10 are given in sec. 10.) The logarithmic term enters the computation 
for n ^ 2 since at least two terms of the series must be used with the above relative error test. Con- 
sequently, using the power series for x ^ 1, this situation will not then be critical when n = \ and 
x ~ 0.56. However, since the code may be used for experimental purposes, the relative error test is 
not applied to the logarithmic term. 

For n ^ 1 and x>l, and n > RMAXI, x > 0, the continued fraction in its "even" form 



E n (x)=e-* 



1 1-71 2(71+1) 



x + n— x + n + 2— x + n + 4— 

= e -*F = e"*EXPENX 

is evaluated in the forward direction. The first convergent F\/Gi = A l /Bi where Ai = l, Am = 
-(M-l)(n + M-2), B M = x + n + 2(M-l). If we define 

F_i = l,Fo = 0, G-j = and G =l 

then successive convergents F mJGm for M= 1, 2, . . . may be obtained by the following recurrence 
relation 

F m = B m F m -\-\-A mF m - 2 

Gm — BmGm -f\-A mGm - 2 . 

Since in the basic continued fraction Am and Bm are always positive with Am^ 1, Fm and Gm 
will always be positive and increasing with increasing M. Care must then be taken to ensure that 
overflow does not occur in generating the successive convergents. For the "even" form of the con- 
tinued fraction, since Am < for M > 1, it is sufficient to check that Gmis always less than RINF/#a/. 
If Gm is greater, thenFM, Fm-i, Gm andG M -i are all scaled by dividing byi?M. 

Since, for the above continued fraction, the successive convergents form a monotonically 
increasing sequence, if through round-off errors, 

l (F M _ JGm - i ) (= PREV) 

1 r , r ^ U \r mi^m ^ r m - lllsM - i ) 

v r mi^m 

then Fm - \\G M - i = F. 

If < 1 - Fm f x ^ n M ~ * ^ TOLER (= 2 " NBM ) , then F m\G m = F. 
r m\Gtm 

The following table gives an indication of the number of terms needed to obtain maximum 
machine accuracy for particular values of NBM, x and n with the various methods of computation. 
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Number 


>f Terms 




Method 


NBM 


= 27 


NBM 


= 60 


x= 1 


n = 1 


n = \0 


n= 1 


n= 10 


Power Series 


13 


14 


21 


22 


Continued Fraction 


26 


17 


113 


76 


Numerical Integration 


129 


513 


16385 


65537 


(Trapezoidal or Simpson's 










Rule) 










x = 22 










Asymptotic Expansion 
Continued Fraction 


22 








5 


6 


12 


16 


Numerical Integration 


2049 








* = 44 










Asymptotic Expansion 
Continued Fraction 


10 


35 


45 




4 


5 


9 


11 


%=70 










Asymptotic Expansion 


8 


15 


23 


62 


Continued Fraction 


4 


5 


8 


10 



4, Range 

For the function e x E „(x), the sum of x and n must be less than or equal to the maximum 
machine value. If the sum is greater, both e x E n {x) and E ,,(x) are set equal to zero. 

For the function E,,(x) the range of x is dominant and essentially equivalent to the range for 
the exponential function. In single precision on the Univac 1108, E „(x) = beyond x ~ 85; in double 
precision beyond x ~ 704. 

5. Accuracy 

Using the Univac 1108 in computing e x E„(x), we find the maximum relative error is 1.3(— 7) 
for the single precision computation and 2.4(— 17) for the double precision computation. 

In computing E n (x), largely due to the error of the exponential routine, the maximum relative 
error is 4(— 7) for the single precision computation and 4.5(— 16) for the double precision 
computation. 

The number of accurate binary digits is essentially the lesser of NBM— NBM 1/2 or NBM— I 
where I is the number of binary digits representing the integer part of x. 

6. Precision 

The precision may be set lower than the maximum by varying NBM or if desirable, deleting 
NBM and setting the proper value of TOLER. 
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7. Timing (Seconds- Univac 1108) 



For n=l, 2, 20 
over the range 


NBM = 


27 


NBM = 60 








* = 0(0.02)1 (153 values) 


0.18 




0.46 


* = 2(1)85 (252 values) 


.19 




.61 


^ = 85(10)715(192 values) 


- 




.28 


Maximum Time/Evaluation 


.003 




.015 



8. Testing 

The double precision results obtained were compared against available published values. 
Further check values were obtained by utilizing multi-precision packages. In addition, double pre- 
cision check values were obtained, where appropriate, with the asymptotic expansion, numerical 
integration, and the "odd" form of the continued fraction. Single precision results were then 
checked against the double precision results. In all cases, the results obtained agreed within the 
reported accuracy. 

9. Driver Program and Its Results 

In the appendix we have included a driver program and its results. The usage of the subroutine 
is thus illustrated and a reasonable set of check values given to simplify the checkout of modifica- 
tions to the subroutine. Two tables are given; one for the function E n (x) and the other for e x E n (x). 
The functions are tabulated to 18 significant figures for n = l, 2, and 20, for x = 0, 10 J (10 J )10 J + 1 
with J = JBEGIN(1) JEND. The value of JBEGIN has been set at -2; the value of JEND is the mini- 
mum value of J for which E\ (x) = 0. 

10. Special Constants 

E t (l) =0.21938 39343 95520 27367 71637 75460 12164 

J (_!)« + !/„. 7l ! = £ 1 (i)+ y = 0.7%59 95992 97053 13428 36758 65542 52408 

i 

^(l)=-y(Euler's constant) 

= -0.57721 56649 01532 86060 65120 90082 40243 
^(2) = l-y 

= 0.42278 43350 98467 13939 34879 09917 59756 
x = e m) = 0.56145 94835 66885 16982 41432 14790 88078 

e* (2) = 1.52620 51115 95863 88047 48887 15036 77561 

e* (3) = 2.51628 68309 39363 58025 62371 75591 18353 

e*<4) = 3.51176 11663 39476 18136 68183 63142 18018 

e* (5) = 4.50919 05949 16874 93725 13380 93818 19453 

e* iG) = 5.50753 78297 01368 13780 61569 70031 98287 

e* {7) = 6.50638 71643 69172 08390 16850 25480 57156 

e*<8) = 7.50554 04760 51118 60659 92002 80578 56992 

e*0) = 8.50489 15798 67776 22381 29291 76413 19981 

e*<io) = 9.50437 85180 84354 74744 37262 91500 55559 

log e 2 =0.69314 71805 59945 30941 72321 21458 17656 
logelO =2.30258 50929 94045 68401 79914 54684 36420 
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2-24 = 0.59604 64477 53906 25 (-7) 
2- 27 = .74505 80596 92382 8125 (-8) 
2~ 36 = .14551 91522 83668 51806 64062 5 (-10) 
2-48= .35527 13678 80050 09293 55621 33789 0625 (-14) 
2-60= 86736 1737 9 88 403 54720 59622 40695 95336 (-18) 
2-io8= .30814 87911 01957 73648 89564 70813 58837 (-32) 

Maximum and Minimum Machine Values and Their Natural Logarithms 
NBC = Number of binary digits in the (biased) characteristic of a floating point number 

2_( 2 NBC-l + i) ^ x < 2 2NBC-1_! 

NBC =8 

2 127 = 0.17014 11834 60469 23173 16873 03715 88410 (39) 
2-129 = 0.14693 67938 52785 93849 60920 67152 78070 (-38) 
log,(2 127 ) =88.02969 19311 13054 29598 84794 25188 42414 
logr(2 129 ) --89.41598 62922 32944 91482 29436 68104 77728 

NBC = 11 

2io23 = o.89884 65674 31157 95386 46525 95394 51236 (308) 
2 M)25 = .27813 42323 13400 17288 62790 89666 55050 (-308) 
log e (2 1023 ) = 709.08956 57128 24051 53382 84602 51714 62914 
log P (2- 1025 ) =-710.47586 00739 43942 15266 29244 94630 98227 
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Appendix 
Implementing Program 



1 
2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 
31 
32 
33 
34 
35 
36 
37 
38 
39 
40 
41 
42 
43 
44 
45 
46 
47 
48 
49 
50 
51 
52 
53 
54 
55 



C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 



SUBROUTINE EXPINT (RN»X»ENXrE'XPENX» IERR) 

LANGUAGE. AMERICAN NATIONAL STANDARD FORTRAN 

DEFINITIONS* 

EN(X)= INTEGRAL (EXP (-X*T) DT/(T**N) ) t FROM 1 
RN(=N)> POSITIVE INTEGER 
X* REAL AND POSITIVE 

SPECIAL CASES 
ENCO)=INFINITY(=RINF* MAXIMUM MACHINE VALUE) 
EN(0)= 1/N-l N .GT. 1 



TO INFINITY 



N .LE« 



E0(X)= EXP(-X)/X 
E0(X)= INFINITY 



X .GT. 1/RINF 
X .LE.' 1/RINF 



USAGE. 



CALL EXPINT (RN> XrENX fEXPENXf IERR) 



(SAME TYPE AS RN) 
(SAME TYPE AS X) 
(SAME TYPE AS X) 



INPUT 

INPUT 
OUTPUT 
OUTPUT 
OUTPUT 



FORMAL PARAMETERS 
RN REAL OR DOUBLE PRECISION TYPF 

FOR A POSITIVE INTEGER N 
X 

ENX= EN(X) 

EXPENX= EXP(X)*EN(X) 
IERR INTEGER VARIABLE 

NORMAL RETURN IERR= 

ERROR RETURN IERR= If X AND/OR N NEGATIVF 

ENX=EXPENX=-INFINITY (IMPOSSIBLE VALUE) 

IERR= 2» N NON-INTEGER 
ENX=EXPFNX=INFINITY 

MODIFICATIONS 
DOUBLE PRECISION UNIVAC 1108 RESULTS ARE OBTAINED IF AS 
SET UP BELOW WHERE 

NBM=ACCURACY DESIRED OR MAXIMUM NUMBER OF BINARY 

DIGITS IN THE MANTISSA 0^ A FLOATING POINT NUMBER 

WITH 

(1) THE DOUBLE PRECISION TYPE STATEMENT 

(2) THE MAXIMUM MACHINE VALUE AMD THE MAXIMUM INTEGER 
<=RMAXI) CONVERTIBLE TO A FLOATING POINT NUMBER 
GIVEN AS DOUBLE PRECISION CONSTANTS 

(3) DOUBLE PRECISION DECIMAL CONSTANTS 

(4) DATA STATEMENT NRM=60 POR THE CONTROL VARIAPLE 

(5) FUNCTION TYPE STATEMENTS - DEXPr DLOG. 

SINGLE PRECISION UNIVAC HOB RESULTS ARE OBTAINED BY 

(1) DELETING THE DOUBLE PRECISION TYPF STATEMENT 

(2) ADJUSTING MAXIMUM MACHINE VALUE 

(3) CHANGING THE D'S TO E f S ON THE DATA CARDS FOR ALL 
DECIMAL CONSTANTS 

(4) SETTING NBM=27 IN THE CONTROL VARIABLE DATA 

(5) CHANGING FUNCTION TYPE - EXP* LOG. 
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56 C FOR. OTHER COMPUTERS THE APPROPRIATE VALUE OF MPM MUST 

57 C BE INSERTED AMD ALL VALUES ADJUSTED ACCORDINGLY, 

58 C 

59 C IF A PRECOMPILED VALUE OF TOLER (=2** (-NBM) ) IS 

60 C INCLUDED IN A DATA STATEMENT* COMPUTATION OF THE CONTROL 

61 C VARIABLE MAY BE OMITTED AND THE RATA STATEMENT FOR NBM 

62 C DELETED. 

63 C 

64 C CAUTION - THE SUBROUTINE CANNOT PEADILY RE ARAPTED TO 

65 C COMPUTE THE EXPONENTIAL INTEGRAL FOR A COMPLEX 

66 C ARGUMENT AS THE CONTINUED FRACTION IS INVALID 

67 C ALONG THE NEGATIVE REAL AXIS. TN ADDITTON MANY 

68 C OF THE COMPARISONS BECOME MEANINGLFSS. 

69 C 

70 C METHOD. 

71 C POWER SERIES* X .LE. K=ULPSr UPPFP LIMIT FOR POWER 

72 C SERIES) 

73 C ENX = SUM(TM)r M=0 1 1 * 2» . . . 'K 

74 C TM =-( (-X)**M/1*2*3...*0/(m-N+1=D) M .ME. M— 1 

75 C = PTERM/D 

76 C TM = PTERM*(LOG(X)-PSI(N) ) M .EO. N-l 

77 C PTERM(0)=-1 

78 C PTERM(M+1>= PTERM ( M) * ( -X ) / ( M+l ) 

79 C PSI(N)= -EULER+1+1/2+ ... +1/(N-1) 

80 C IP R.E. (=ARS(TM/SUM) ) .LE. TOLER ( =?** (-NRM) ) , V!=K 

81 C 

82 C CONTINUED FRACTION* X .GT. 1 

83 C EN(X)=EXP(-X)*(1I/I(X+N)- 1 *N I/I ( X+N+?) - 

84 C 2*(M+l)l/i(X+N+<0-... 

85 C EN(X)=EXP(-X)*II(AM I/I BM> M=1,2,...,K 

86 C AM(1)=1 AM(M)=-(M-D*(N+M-2) 

87 C BM(M)=x+N+2*(^-l) 

88 C 

89 C EN(X)=EXP(-X)*FM(K)/GM(K)=EXP(-X)*F(K) 

90 C IF R.E. UABSQ-PRFV/F) ) .LE. TOLEP (=2** (-NRM) ) t M=K 

91 C F=FM/GM PREV=FMM1/GMM1 

92 C 

93 C RANGE. 

94 C FOR EXP(X)*EN(X) (N+X) .LE. MAXIMUM MACHINF VALUE 

95 C FOR EN(X) X=APPROXIMATELY *5.0r SINGLE PRECISION 

96 C 704.0* DOUBLE PRECISION 

97 C BEYOND THIS RANGE EN(X>=0 

98 C 

99 C ACCURACY. THE NUMBER OF ACCURATE BINARY DIGITS IS ESSEN- 

100 C TIALLY THE LESSER OF 

101 C NBW - SQUARE ROOT OF NBM 

102 C NBM - I (NUMBER OF BINARY DIGITS REPRESENT- 

103 C IMG THE INTEGER PART OF X) 

104 C (ACCURACY OF THE EXPONENTIAL ROUTINE) 

105 C 

106 C PRECISION. VARIABLE - BY SETTING THE DESIRED MRM OR TOLER. 

107 

108 C MAXIMUM TIME. UNIVAC 1108-EXEC II 

109 C (SECONDS) NRV|=27 NRM=60 

110 C .0025 .015 

111 C 

112 C STORAGE. COMPILED BY THE UNIVAC llOfW EXEC R/FORTRAN Vr 

113 C THIS SUBROUTINE REQUIRES 478 WORDS OF STORAGE. 

114 C 

115 C 

116 C 
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117 C 

118 C 

119 C 

120 C 

121 C 

122 C MACHINE DEPENDENT STATEMENTS 

123 C 

124 DOUBLE PRECISION AM* ASUMr BMt D* FNX* EULER* 

125 1 EXPENX f EXPNX * F * FM * FMM1 * FMM2 * GM * GMM1 , 

126 2 GMM2 * HALF r ONE f ONPT^V * PREV r PS I * PTERM r 

127 3 RINFfRMfRMAXIfRNfRNMlrPRNfSUMrTEMPt 

128 4 TM * TOLER* TWO *ULPS*X*XLOG*ZFRO 

129 C 

130 C CONSTANTS 

131 C MAXIMUM MACHINE VALUE 

132 C MAXIMUM CONVERTIBLE INTFGER 

133 C 

134 DATA RINF / .898846567431157954D308 / 

135 DATA RMAXI / 134217727.D0 / 

136 C 

137 C ADDITIONAL CONSTANTS 

138 C 

139 DATA EULER / .577215664901532R61D0 / 

140 DATA HALF p ONE fONPTFVr TWO tULPSr ZERO / 

141 1 •5D0*1.D0*1.5D0*2.D0*1,D0*0.0D0 / 

142 C 

143 C CONTROL VARIABLE 

144 C 

145 C NBM=ACCURACY DESIRED OR MAXIMUM NUMBER OF BINARY 

146 C DIGITS IN THE MANTISSA OF A FLOATING POINT NUMBER 

147 C 

148 DATA NBM / 60 / 

149 TOLER=TWO**(-NBM) 

150 C 

151 C VALIDITY TEST FOR INPUT PARAMETERS 

152 C ERROR CONDITION* IMPOSSIBLE VALUES RFTURNED 

153 C 

154 C NEGATIVE ZERO CHECKS 

155 C 

156 IERR=0 

157 IF (RN .LT. ZERO) GO TO 3 

158 IF (RN .GT. RMAXI) GO TO 10 

159 C 

160 C VALIDITY TEST FOR N INTEGER 

161 C 

162 N=RN 

163 RRN=N 

164 IF (RRN) 2*1*2 

165 1 IF ((RN-RRN)-TOLER) 10*10*2 

166 2 IF (<ONE-(RRN/RN>> .LE* TOLER) GO TO 10 

167 RRN=RRN+ONE 

168 IF ( (ONE-(RN/RRN)) .LE. TOLER) GO TO 10 

169 C 

170 C ERROR RETURN 

171 C N NON-INTEGER 

172 C 

173 IERR=2 

174 ENX=RINF 

175 EXPENX=ENX 

176 RETURN 

177 C 
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178 




3 


IF (-RN) 12*10*12 


179 








180 




10 


IF (X) 11*30*20 


181 




11 


IF (-X) 12*30*12 


162 


c 






183 


c 




ERROR RETURN 


184 


c 




X AND/OR N NEGATIVF 


185 


c 






186 




12 


IERR=1 


187 






ENX=-RINF 


188 






EXPENX=ENX 


189 






RETURN 


190 


c 






191 


c 




FUNCTION TYPE STATEMENTS 


192 


c 






193 




20 


EXPNX=DEXP(-X) 


194 






XLOG=OLOG(X) 


195 


c 






196 






IF (RN .GE. HALF) GO TO 60 


197 


c 






198 


c 




SPECIAL CASES 


199 


c 






200 






IF (X .LE. ONE/RINF) GO TO 40 


201 






EXPENX=ONE/X 


202 






ENX=EXPNX*EXPENX 


203 






RETURN 


204 


c 






205 




30 


1^ (RN .LT. ONPTFV) GO TO 40 


206 






ENX=ONE/(RN-ONE> 


207 






GO TO 50 


208 


c 






209 




40 


ENX=RINF 


210 




50 


EXPENX=ENX 


211 






RETURN 


212 


c 






213 




60 


IF (RN .GT. RMAXI) GO TO 70 


214 






IF (X .LE. ULPS) GO TO 100 


215 




70 


IF (X .LE. (RINF-RN)) GO TO 200 


216 






EXPENX=ZERO 


217 






ENX=EXPENX 


218 






RETURN 


219 


c 






220 


c 




METHOD POWER SERIES 


221 


c 




v 


222 


100 


RM=0 


223 






PTERM=-ONE 


224 






SUM=0 


225 






PSI=-EULER 


226 






D=-(RN-ONE> 


227 


c 






228 


110 


IF (D .GE. HALF) GO TO 130 


229 






IF (-D .GE. HALF) GO TO 120 


230 


c 






231 


c 




COMPUTE TM FOR m .EQ. N- 


232 


c 






233 






SUM=PT£RM* ( XLOG-PSI ) +SUM 


234 






GO TO 170 


235 


c 






236 


c 




COMPUTE PSI (N) 


237 


c 






238 


120 


PSI=PSI+(ONE/(RM+ONF) ) 



209 



239 


C 




240 


c 


COMPUTE TM FOR M .NE. N-l 


241 


c 




242 


130 


TM=PTERM/D 


243 




SUM=TM+SUM 


244 


C 




245 


C 


TOLERANCE CHECK 


246 


C 




247 


C 


ZERO CHECKS 


248 


c 




249 




IF (SUM .LT. ZERO) GO TO 140 


250 




ASUM=SUM 


251 




GO TO 150 


252 


c 




253 


140 


ASUMr-SUM 


254 


150 


IF (ASUM) 160*170*160 


255 


160 


IF (TM .LT. ZERO) TM=-TM 


256 




IF ( TM/ASUM .GT. TOLER) GO TO 


257 


C 




258 




ENX= SUM 


259 




EXPENX=ENX/EXPNX 


260 




RETURN 


261 


C 




262 


C 


ADDITIONAL TERMS 


263 


C 




264 


170 


RM=RM+ONE 


265 




D=D+ONE 


266 




PTERM=-(X*PTERM)/RM 


267 




GO TO 110 


268 


c 




269 


c 


METHOD -*- CONTINUED FRACTION 


270 


c 




271 


200 


RM=ONE 


272 




FMM2=0NE 


273 




GMM2=0 


274 




FMM1=0 


275 




GMM1=0NE 


276 




PREV=FMM1/GMM1 


277 




AM=0NE 


278 




BM=X+RN 


279 




RNM1=RN-0NE 


280 


C 




281 


210 


FM=BM*FMM1 + AM+FMM2 


282 




GM=BM*GMM1 + AM*GMM2 


283 




F=FM/GM 


284 


C 




285 


C 


TOLERANCE CHECK. 


286 


C 




287 




TEMP=ONE-(PREV/F) 


288 




IF (TEMP •GT. ZERO) GO TO 220 


289 


C 




290 




EXPENX=PREV 


291 




GO TO 230 


292 


C 




293 


220 


IF (TEMP .GT. TOLER) GO TO 240 


294 


C 




295 




EXPENX=F 


296 


230 


ENX=EXPNX*EXPENX 


297 




RETURN 


298 


C 




299 


240 


IF( GM .LT. RINF/BM) GO TO 250 



210 



300 C 

301 C SCALING 

302 C 

303 C BOTH FM AND GM MUST BE TESTED IF N=l AND X .LT. .11 

301 C SCALING SHOULD NOT BE DELETED AS THE VALUES OF FM AND GM 

305 C MAY OVERFLOW FOR SOME CHOICES OF THE PARAMFTERS. 

306 C 

307 FMM1=FMM1/BM 

308 GMM1=GMM1/BM 

309 FM=FM/BM 

310 GM=GM/BM 

311 C 

312 C ADDITIONAL CONVERGENTS 

313 C 

311 250 AM=-RM*(RNM1+RM) 

315 RM=RM+ONE 

316 BM=BM+TWO 

317 FMM2=FMM1 

318 GMM2=GMM1 

319 FMM1=FM 

320 GMM1=GM 

321 PREV=F 

322 GO TO 210 

323 C 

321 END 



211 



Driver (Test) Program 

1 C 

2 C DOUBLE PRECISION i;EST PROGRAM 

3 C ALGORITHM FOR EXPONENTIAL INTEGRAL EN(X) AND EXP(X)*EN(X) 

4 C RN(=NrA POSITIVE INTEGER)* Xr REAL AND POSITIVE 

5 C (AMERICAN NATIONAL STANDARD FORTRAN - UNIVAC 1108) 

6 C 

7 C USAGE OF CALL EXPINT (RNrXtENXt EXPENX* IERR) 

8 C 

9 C TABLE OF ENX (X(L) p N(K) ) AND EXPENX <X(L) rN(K) ) 

10 C JT=1 JT=2 

11 C L=1»2*.».>NCX* K=1*2*...*NN 

12 C N=l*2*20 X=0» 10**J (10**J) 10**(J+1) J=JBEGIN( 1 ) JEND 

13 C (10**J) (ID=1*2* ...*9) JBEGIN=-2 

14 C JEND=MIN J FOR ENX <X= ( 10**J) *ID#N( 1 ) )=0 

15 C 

16 C NLINE = NUMBER OF AVAILABLE PRINT" LINES PER PAGE 

17 C NFCOL = NUMBER OF FUNCTIONAL COLUMNS PER PAGE 

18 C LINHD = NUMBER OF HEADING LINES PER PAGE 

19 C 

20 C RMAXI = MAXIMUM CONVERTIBLE INTEGER 

21 C 

22 C DIMENSION N(NN) #X(NX) * ENX(NXt NN) » EXPENX ( NX »NN) 

23 C 

24 C MACHINE DEPENDENT STATEMENTS 

25 C INPUT DATA 

26 C 

27 DOUBLE PRECISION D>OFCTR* DX , ENX* ETPTFW EXPENXf 

28 1 PTFIVEfRMAXIrRNrTENfXrZERO 

29 C 

30 DIMENSION N(3) rX(100) * ENX ( 100 p 3) p FXPENX C 100 p 3) 

31 C 

32 DATA RMAXI / 134217727. UO / 

33 DATA ETPTFV,PTFIVE*ZERO / 8.5D0 p . 5D0r 0. 0D0 / 

34 C 

35 DATA NN*NX>NLINE>NFCOL*LINHD/ 3*100*58*3*3 / 

36 DATA N(l) *N(2) *N(3) / 1*2*20 / 

37 DATA JBEGIM / -2 / 

38 C 

39 1000 FORMAT < 1H1 p 38X p 5HEN( X ) // ) 

40 1001 FORMAT (4X* 1HX* 9X *2HN=* II 12 *4 ( 13X* I 12) ) 

41 1002 FORMAT ( 1D7. 1 * 5025. 18) 

42 1003 FORMAT ( 1H1 p 35X p 12HEXP (X) *EN (X) //) 

43 C 

44 C FUNCTION TYPE STATEMENT 

45 C 

46 MAXJI=DLOG10 (RMAXI ) 

47 C 

48 C COMPUTATION AMD STORAGE OF FUNCTIONS 

49 C 

50 C SET UP INITIAL X 

51 C 

52 D=ZERO 

53 C 

54 C SET UP INITIAL DECADE 

55 C 

56 TEN=10 

57 L=l 

58 J=JBEGIN 

59 ID=1 
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! 


60 


c 








61 






IF (IABS(J) .GT. MAXJT) GO TO 1 




62 






DFCTRzlO**(IABS(J)) 




63 






GO TO 2 




64 


c 








65 




1 


DFCTR=TEN**(IABS(J>) 




66 


c 








67 




2 


IF (J .LE. 0) GO TO 3 




68 






DX=D*DFCTR 




69 






GO TO 4 




70 


c 








71 




3 


DX=D/DFCTR 




72 


c 








73 




4 


X(L)=DX 




74 


c 








75 






DO 5 K=lrNN 




76 






RN=N(K) 




77 




5 


CALL EXPINT(RNfOXtFNX(LfK) >EXPFNX(L*K) * IERR) 




78 


c 








79 






IF (ENX(Lrl) .LE. ZERO) GO TO 9 




80 






IF (L .GE. NX) GO TO 9 




81 






IF (D .LT. PTFIVE) GO TO 8 




82 






IF (D .GT. ETPTFV) GO TO 6 




83 






ID=ID+1 




84 






GO TO *8 




85 


c 








86 




6 


J=J+1 




87 






ID=1 




88 






IF (J .LE. 0) GO TO 7 




89 






DFCTR=DFCTR*TEN 




90 






GO TO 8 




91 


c 








92 




7 


DFCTR=DFCTR/TEN 




93 


c 








94 




8 


D=ID 




95 






L=L+1 




96 






GO TO 2 




97 


c 








98 




9 


NCX=L 




99 


c 








100 


c 




PRINTING OF TABLES 




101 


c 




L=lr2f ...fNCX 




102 


c 




K=lf2r...fNN 




103 


c 








104 






JT=1 




105 


c 








106 




21 


NA=1 




107 






NB=NFCOL 




108 


c 








109 




22 


IF (NB .GT. NN) NB=NN 




110 






L=l 




111 






IF (X(l)) 24*23*24 




112 




23 


NPLINE=1 




113 






LCOUNT=NLINE-LINHD-l 




114 






INDEXS=9 




115 






GO TO 26 




116 


c 








117 




24 


NPLINE=0 




118 


c 








119 




25 


LCOUNT=NLINE-LINHD 




120 






INDEXS=1 
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121 26 NPLINE=NPLINE+9*(LCOUNT/10) 

122 IF (JT .EQ. 1) GO TO 27 

123 WRITE (6*1003) 

124 GO TO 28 

125 C 

126 27 WRITE (6*1000) 

127 28 WRITE (6*1001) (N(K) *K=NA* NB) 

128 C 

129 29 IF (JT .EQ. 1) GO TO 30 

130 WRITE (6*1002) X(D * (EXPENX(L*K) * K=NA*NB) 

131 GO TO 31 

132 C 

133 30 WRITE (6*1002) X(L) * (FNX(L*K) *K=NA*NB) 

134 31 IF (L .GE. NCX) GO TO 35 

135 IF (L .GE. NPLINE) GO TO 34 

136 IF (INDEXS .EQ. 9) GO TO 32 

137 INDEXS=INDEXS+1 

138 GO TO 33 

139 C 

140 32 WRITE (6*1004) 

141 INDEXS=1 

142 C 

143 33 L=L+1 

144 GO TO 29 

145 C 

146 C EXCESS X*S 

147 C 

148 34 L=L+1 

149 GO TO 25 

150 C 

151 C EXCESS N'S 

152 C 

153 35 IF (NB .GE. NN) GO TO 36 

154 NA=NA+NFCOL 

155 NB=NB+MFCOL 

156 GO TO 22 

157 C 

158 C ADDITIONAL FUNCTION 

159 C 

160 36 IF (JT .GE. 2) STOP 

161 JT=JT+1 

162 GO TO 21 

163 C 

164 1004 FORMAT (1H ) 

165 C 

166 END 
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Test Results 



EN(X) 



X N= 1 

.0 .898846567431157951+308 

.1-001 .40379295765381 1384+001 

.2-001 •335470778330970952+001 

.3-001 .29591 IB 72402 128069+001 

.4-001 .268126368902527992+001 

.5-001 .246789848850997437+001 

.6-001 .229530691814378233+001 

.7-001 .215083818025679885+001 

.8-001 .202694100258574148+001 

.9-001 .191874477003266286+001 

.1+000 .182292395841939067+001 

.2+000 .1222650 54418389309+001 

.3+000 .905676651675846714+000 

.4+000 .702380118865662480+000 

.5+000 .559773594776160810+000 

.6+000 .454379503189402111+000 

.7+000 .373768843233509144+000 

.8+000 .310596578545543035+00 

.9+000 .260183939325999640+000 

.1+001 .219383934395520273+000 

.2+001 .489005107080611205-001 

.3+001 .130483810^41970376-001 

.4+001 .377935240984890646-002 

.5+001 .114829559127532581-002 

.6+001 .360082452162658656-003 

.7+001 .115481731610338216-003 

.8+001 .376656228439249018-004 

.9+001 .124473541780062721-004 

.1+002 .415696892968532427-005 

.2+002 .983552529064988168-010 

.3+002 .302155201068881253-014 

.4+002 .103677326145165697-018 

.5+002 .378326402955045901-023 

.6+002 .143586756568125679-027 

.7+002 .560030628581343844-032 

.8+002 .222854325868847290-036 

.9+002 .9005474 05886741103-041 

.1+003 .368359776168203217-045 

.2+003 .688522610630763559-089 

.3+003 .171038427680451011-132 

.4+003 .477601358642097222-176 

.5+003 .142207678225363842-219 

.6+003 .440998979450983796-263 

.7+003 .140651876623403292-306 

.8+003 .000000000000000000 



2 20 

.1000 00000000 000000 + 001 .526315789473684210-001 



.94967053798378691 6+000 
.913104517640561112+000 
.881671971827869758+000 
.853538891591312014+000 
.827834500075215292+000 
.804046118495621770+000 
.781835147287972310+000 
.760961066179776466+0 00 
.741244155968288530+000 

.722545022194020509+000 
.574200644241203245+000 
.469115225178963854+000 
. 389367998489374312+000 
.32664386232455301.7+000 
.276183934180 385169+00 
.234947113527953114+000 
.200851701280 787165+00 
. 1724041143471°9437+000 

. 148495506775922048+000 
.375342618204904530-001 
.1064192508527P8306-0 01 
.319822924933855437-0 02 
.996469042708838118-0 03 
.318257463690406465-0 03 
.10350984428214R693-0 03 
.341376451511126247-004 
.113836164846231004-004 

.383024046563160875-005 
.940485643085814897-010 
.2Q296693677373696 Q -014 
.101261209484^61107-018 
.371178331886882736-023 
. 141305368608979607-027 
.552353358392398957-0 3? 
.220167808946368433-0 36 
.890859710098454963-041 

.364782143388037826-045 
.685130547521041110-089 
.170473919984834340-132 
.476416214561680312-176 
.141924954730934210-219 
.440267629840803335-263 
.140451801215403972-306 
.000000000000000000 



.520789541793351476-001 
.515321496513528285-001 
.509911038545502872-001 
.5045^7559326025865-001 
.499260456747777298-001 
.494019135090578269-001 
.488833004953339252-001 
.483701483186737100-001 
.478623992826612893-001 

.473599963028082897-001 
.426178656910138496-001 
.383518490868851917-001 
.345140002098048623-001 
.310612173936309823-001 
.279547522194935671-001 
.251597681997253041-001 
.226449443896666456-001 
.203821193314687553-001 

.183459712067558733-001 

.641430585532489937-002 
.2248648064140 00919-002 
.790189099803715192-003 
.278274592885730814-003 
.981884227302490187-004 
.347068486248843707-0 04 
.122877543566663061-004 
.435687424044296172-005 

.154693627987772485-005 
.521648146507963053-010 
. 188625975171560854-014 
.711928232079808894-019 
.276642308097688871-023 
.109793192761257778-027 
.442791440616991968-032 
.180841147365357070-036 
.746125498994224778-041 

.310431603951609625-045 
.629301793252036997-089 
. 160912502553792433-132 
.456044233853820537-176 
.137021182168391364-219 
.427505488621821494-263 
.136945221165125589-306 
.000000000000000000 
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EXP(X)*EN(X) 



X N= 1 

.0 .898846567431157951+308 

.1-001 .407851144345642585+001 

.2-001 .342247737593075321+001 

.3-001 .304923730567447425+001 

.4-001 .279068813598834047+001 

.5-0 01 .259443034976061332+001 

.6-001 .243724077122346628+001 

.7-001 .230679154487934795+001 

.8-001 .219575897504124847+001 

.9-001 .209944118436360746+001 

.1+00 .201464254470845168+001 

.2+00 .149334874693223961+001 

.3+00 .122253560508058556+001 

.4+000 .104782800845600643+001 

.5+00 .922910632483730466+000 

.6+000 .827933435273508820+000 

.7+00 .752678020029587137+000 

.8+00 .691245397802831489+000 

.9+00 .639949226639299739+000 

.1+001 .596347362323194074+000 

.2+001 .361328616888222592+000 

.3+001 .262083740255318501+000 

.4+001 .206345649901055832+000 

.5+001 .170422176284732204+000 

.6+001 .145267629233886893+00 

.7+001 .126641096076632765+000 

.8+001 .112279639253499312+000 

.9+001 .100861955580640929+000 

.1+002 .915633339397880818-001 

.2+002 .477185454959608417-001 

.3+002 .322897387589801252-001 

.4+002 .244041150796285762-001 

.5+002 .196151099301148704-001 

.6+002 .163977137080465268-001 

.7+002 .140872270003268122-001 

.8+002 .123475166636786097-001 

.9+002 .109903102083356454-001 

.1+003 .990194228673301838-002 

.2+003 .497524632317935662-002 

.3+003 .332229556527070706-002 

.4+003 .249378101793988503-002 

.5+003 .199601590476041089-002 

.6+003 .166389810215794723-002 

.7+003 .142653641830088669-002 

.6+003 .124844139167435033-002 



2 20 

.100000000000000000+001 .526315789473684210-001 



.959214885565435741+000 
.931550452481384936+000 
.908522880829765773+000 
. 888372474560466383+000 
.870278482511969335+000 
.853765553726592023+000 
.838524591858445643+000 
.824339281996700124+000 
.811050293407275329+000 

.798535745529154834+000 
.701330250613552081+000 
.633239318475824335+000 
.580868796617597430+000 
.538544683758134765+000 
.503239938835894713+000 
.473125385979289004+000 
.447003681757734811+000 
.424045696024630237+000 

.403652637676805926+000 
.277342766223554833+000 
.213748779234044509+000 
. 174617400395776667+000 
.147889118576338992+000 
. 128394224596678636+000 
.113512327463570647+000 
.101762885972005506+000 
.922423997742316384-001 

.843666606021191810-001 
.456290900807831660-001 
.313078372305962434-001 
.238353968148569492-001 
. 192445034942564817-001 
.161371775172083932-001 
.138941099771231436-001 
.121986669057112264-001 
.108720812497919106-001 

.980577132669815934-002 
.495073536412867515-002 
.331133041878788067-002 
.248759282404598474-002 
.199204761979455498-002 
.166113870523165916-002 
.142450718937931575-002 
. 124688666051973815-002 



.526023563704061981-001 
.525731681287694728-001 
.525440141582159931-001 
.525148943946742833-001 
.524858087742430369-001 
.524567572331905115-001 
.524277397079539275-001 
.523987561351388700-001 
.523698064515186901-001 

.523408905940339148-001 
.520535787019039881-001 
.517695812767576684-001 
.514888379273828789-001 
.512112898107201955-001 
.509368795793975367-001 
.506655513315264279-001 
.503972505626390015-001 
.501319241196527574-001 

.498695201567573519-001 
.473956658006950287-001 
.451653037195386256-001 
.431428630256384291-001 
.412996114281546988-001 
.396120369170190107-001 
.380606810266053526-001 
.366292794923112188-001 
.353041176322005723-001 

. 340735390554722559-001 
.253085524935690803-001 
.201574668908614936-001 
. 167577416876913709-001 
. 143430890423517574-001 
. 125384637484814277-001 
.111381471287998954-001 
.100197250911190282-001 
.910574016947094732-002 

.834476515943799403-002 
.454731825027371299-002 
.312560692284812426-002 
.238122114425943504-002 
. 192321864972550092-002 
. 161298688732910590-002 
. 138894233048444955-002 
. 121954838359794140-002 
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